Identification of Unique Key miRNAs, TFs, and mRNAs in Virulent MTB Infection Macrophages by Network Analysis

Although Mycobacterium tuberculosis (MTB) has existed for thousands of years, its immune escape mechanism remains obscure. Increasing evidence signifies that microRNAs (miRNAs) play pivotal roles in the progression of tuberculosis (TB). RNA sequencing was used to sequence miRNAs in human acute monocytic leukemia cells (THP-1) infected by the virulent MTB-1458 strain and the avirulent vaccine strain Mycobacterium bovis Bacillus Calmette-Guérin (BCG). Sets of differentially expressed miRNAs (DE-miRNAs) between MTB-1458/BCG-infected groups and uninfected groups were identified, among which 18 were differentially expressed only in the MTB-1458-infected THP-1 group. Then, 13 transcription factors (TFs) and 81 target genes of these 18 DE-miRNAs were matched. Gene Ontology classification as well as Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis showed that the candidate targets were predominantly involved in apoptotic-associated and interferon-γ-mediated signaling pathways. A TF-miRNA-mRNA interaction network was constructed to analyze the relationships among these 18 DE-miRNAs and their targets and TFs, as well as display the hub miRNAs, TFs, and target genes. Considering the degrees from network analysis and the reported functions, this study focused on the BHLHE40-miR-378d-BHLHE40 regulation axis and confirmed that BHLHE40 was a target of miR-378d. This cross-talk among DE-miRNAs, mRNAs, and TFs might be an important feature in TB, and the findings merited further study and provided new insights into immune defense and evasion underlying host-pathogen interactions.


Introduction
Tuberculosis (TB) is a zoonotic communicable disease caused by Mycobacterium tuberculosis complex (MTBC) and is still one of the deadliest diseases worldwide. According to the

Expression Profiles of miRNAs in THP-1 Cells after Infection
miRNA expression profiles of MTB-1458-infected THP-1 cells or BCG-infected THP-1 cells at 6 and 24 h postinfection (hpi) were constructed. Uninfected THP-1 was used as the control. Generally, the expression profile of BCG-infected groups was more similar to uninfected groups than MTB-1458-infected groups ( Figure 1A). Clean reads mapping against human mature miRNAs of miRBase results demonstrated that the proportions of miRNA reads in MTB-1458 groups were only 35% (6 hpi) and 33% (24 hpi), which were dramatically lower than in BCG and uninfected groups (all >70%; Table 1). These indicated that MTB-1458 induced a stronger miRNA expression profile variation and inhibited miRNA expression in the entire level. upregulated at 6 hpi, while miR-212-3p and miR-132-3p were both upregulated at 24 hpi. Fourteen miRNAs at 6 hpi and thirteen miRNAs at 24 hpi were only differentially expressed after BCG infection, of which six miRNAs overlapped. Fifty miRNAs at 6 hpi and fifty-five miRNAs at 24 hpi were specifically regulated by MTB-1458, of which thirty-five miRNAs overlapped ( Figure 1C; Table S1). Taken together, the results indicated that DE-miRNAs had little differences between 6 and 24 hpi in either BCG or MTB-1458 groups but showed a distinct difference within virulent strain-infected groups and avirulent strain-infected groups.  A total of 1685 mature miRNAs were identified, including 234 novel miRNAs. For the 1451 known miRNAs, 92 were differentially expressed between MTB-1458-infected cells or BCG-infected cells and uninfected cells with a fold change of ≥2 or ≤0.5 (p ≤ 0.01; Figure S1). In detail, there were 22 DE-miRNAs between BCG-infected cells and uninfected cells, 15 upregulated DE-miRNAs at 6 and 15 upregulated DE-miRNAs at 24 hpi, and 8 miRNAs overlapped. There were 73 DE-miRNAs between MTB-1458-infected and uninfected cells. DE-miRNAs were 51 at 6 hpi, of which 43 miRNAs were downregulated and 8 were upregulated. At 24 hpi, DE-miRNAs were 57, of which 42 were downregulated and 15 were upregulated. Among them, 28 downregulated and 7 upregulated miRNAs were both differentially expressed at 6 and 24 hpi ( Figure 1B,C).
In both the MTB-1458-infected group and the BCG-infected group, miR-146a-3p was upregulated at 6 hpi, while miR-212-3p and miR-132-3p were both upregulated at 24 hpi. Fourteen miRNAs at 6 hpi and thirteen miRNAs at 24 hpi were only differentially expressed after BCG infection, of which six miRNAs overlapped. Fifty miRNAs at 6 hpi and fifty-five miRNAs at 24 hpi were specifically regulated by MTB-1458, of which thirty-five miRNAs overlapped ( Figure 1C; Table S1). Taken together, the results indicated that DE-miRNAs had little differences between 6 and 24 hpi in either BCG or MTB-1458 groups but showed a distinct difference within virulent strain-infected groups and avirulent straininfected groups.

Target and TF Prediction of Unique DE-miRNAs in MTB-1458 Groups
The potential targets of 18 selected DE-miRNAs in MTB-1458-infected groups were predicted by the miRDB and miRanda online databases, together with mRNAseq data. The potential TFs of the selected miRNAs were predicted by the TransmiR database. Fourteen potential targets and twelve candidate TFs were predicted for five upregulated DE-miRNAs, and sixty-seven potential targets and eleven candidate TFs were predicted for thirteen downregulated DE-miRNAs (Table 2). Excluding repeat TFs, 81 potential targets and 13 TFs were selected for further analysis.

Gene Ontology (GO) Classification and Kyoto Encyclopedia of Genes and Genomes (KEGG Pathways) Pathway Enrichment Analysis of the Candidate Targets
To further understand the immune regulatory functions of the specific DE-miRNAs, the final candidate targets were subjected to GO classification and KEGG pathway enrichment analysis. The candidate targets were enriched in iron and protein, as well as DNA-binding-related molecular functions (MF) and interferon-γ (IFN-γ)-mediated signaling pathway and apoptotic-process-associated biological processes (BP), located in the cytosol of cellular component (CC; Figure 3A). For KEGG pathway enrichment analysis, 5 of the top 10 enriched pathways were associated with TB-related immune response, such as the Foxo signaling pathway (target genes: PRKAG2, PLK3, TNFSF10, FBXO32, and SGK1), tumor necrosis factor (TNF) signaling pathway (TRAF1, IRF1, CSF1, and CEBPB), phosphatidylinositol 3-kinase (PI3K-)/Akt signaling pathway (PIK3R5, MYB, NGFR, SGK1, and CSF1), Rap1 signaling pathway (CTNNB1, CSF1, MRAS, and NGFR), and C-type lectin receptor signaling pathway (IRF1, PLK3, and MRAS; Figure 3B). These candidate targets of DE-miRNAs might be pivotal mediators during MTB infection.

Target and TF Prediction of Unique DE-miRNAs in MTB-1458 Groups
The potential targets of 18 selected DE-miRNAs in MTB-1458-infected groups were predicted by the miRDB and miRanda online databases, together with mRNAseq data. The potential TFs of the selected miRNAs were predicted by the TransmiR database. Fourteen potential targets and twelve candidate TFs were predicted for five upregulated DE-miRNAs, and sixty-seven potential targets and eleven candidate TFs were predicted for thirteen downregulated DE-miRNAs (Table 2). Excluding repeat TFs, 81 potential targets and 13 TFs were selected for further analysis.

Gene Ontology (GO) Classification and Kyoto Encyclopedia of Genes and Genomes (KEGG Pathways) Pathway Enrichment Analysis of the Candidate Targets
To further understand the immune regulatory functions of the specific DE-miRNAs, the final candidate targets were subjected to GO classification and KEGG pathway enrichment analysis. The candidate targets were enriched in iron and protein, as well as DNAbinding-related molecular functions (MF) and interferon-γ (IFN-γ)-mediated signaling pathway and apoptotic-process-associated biological processes (BP), located in the cytosol of cellular component (CC; Figure 3A). For KEGG pathway enrichment analysis, 5 of the top 10 enriched pathways were associated with TB-related immune response, such as the Foxo signaling pathway (target genes: PRKAG2, PLK3, TNFSF10, FBXO32, and SGK1), tumor necrosis factor (TNF) signaling pathway (TRAF1, IRF1, CSF1, and CEBPB), phosphatidylinositol 3-kinase (PI3K-)/Akt signaling pathway (PIK3R5, MYB, NGFR, SGK1, and CSF1), Rap1 signaling pathway (CTNNB1, CSF1, MRAS, and NGFR), and C-type lectin receptor signaling pathway (IRF1, PLK3, and MRAS; Figure 3B). These candidate targets of DE-miRNAs might be pivotal mediators during MTB infection.

TF-miRNA-mRNA Network
Combined with the above analysis, the complex regulatory network of miRNAs, mRNAs, and TFs was constructed and analyzed for key miRNAs, mRNAs, and TFs. In the five upregulated DE-miRNAs networks, miR-335-3p had the highest node degree of interactions, with seven targets and six TFs, followed by miR-1290 with four targets and eight TFs, miR-146-5p with zero targets and nine TFs, and miR-96-5p with five targets and four TFs. miR-625-3p had the minimum regulation degree, with one target and six TFs. Three (MYB, SH3KBP1, and TSC22D1) of fourteen genes were targeted by two miRNAs, and others only had one associated miRNA. In TFs, JUND was associated with all five miRNAs, followed by CEBPB and JUN, each targeting four miRNAs. In all interactions, JUND-miR-335-3p-MYB/SH3KBP1 was regarded as the key axis because of the highest degree of all three aspects in the network, which might affect the progress of MTB infection ( Figure 4A).

Key DE-miRNA Validation by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Six DE-miRNAs (logFC > 2 at 6 and 24 hpi) were further validated by qRT-PCR according to the TF-miRNA-mRNA network analysis. miR-3184-5p, miR-320b, miR-140-3p, miR-107, and miR-378d were downregulated, whereas miR-1290 was upregulated in MTB-1458-infected groups compared to uninfected groups, in accordance with sequencing data ( Figure 5A,B). In 13 downregulated DE-miRNAs network, miR-3184-5p had the highest degree of interactions, with 12 targets and 8 TFs, followed by miR-320b with 17 targets and 1 TF, miR-107 with 8 targets and 8 TFs, and miR-140-3p with 9 targets and 7 TFs. Others had <16 connections. miR-378c and miR-4286 had only targets but no TFs predicted because their relative data were absent in the database. In the 67 targets, ZCCHC2 had the highest degree and was regulated by 4 miRNAs, followed by TNPO1, BHLHE40, AFF1, ARL8B, and RDX, all of them are targeted by 3 miRNAs. Eleven out of sixty-seven were regulated by two miRNAs, and fifty of them were targeted by only one miRNA. For TFs, the top ranking was MYB, which was linked to seven miRNAs. BHLHE40, CEBPB, IRF1, RUNX3, and JUND were correlated with five miRNAs. PML, ATF3, and JUN were involved with four miRNAs. In addition, TRIM25 and HNF4G were associated with three and two miRNAs, respectively ( Figure 4B).
Considering the degrees from network analysis and reported functions, this study found that (i) MYB-miR-3184-5p-BMF/IRF1/PML had the highest degree of all TFs, miRNA, and mRNAs in the network. (ii) Three miRNA interactions, (IRF1-miR-3184-5p-IRF1, CEBPB-miR-3184-5p-CEBPB, and BHLHE40-miR-107/miR-378d-BHLHE40), presented the same targets and TFs, respectively, and the whole connection formed like a loop. (iii) Some targets, such as ARL8B, and TFs, such as MYB, IRF1, and CEBPB, participated in the TB immune-related BPs and had a higher degree in the network. However, these findings need to be validated in a series of subsequent experiments ( Figure 4B).
qRT-PCR results indicated that BHLHE40 mRNA expression was significantly increased in MTB-1458-infected THP-1 macrophages in comparison to uninfected cells and was negatively correlated with the miR-378d expression trend ( Figure 6A). The interaction between miR-378d and BHLHE40 was identified by the dual-luciferase reporter using psiCHECK-BHLHE40-wild and psiCHECK-BHLHE40-mut plasmids ( Figure 6B). Results showed the Renilla/firefly ratio was significantly lower in miR-378d mimic than in negative control (NC) mimic in the psiCHECK2-BHLHE40-wild group, and the inhibition ef-
qRT-PCR results indicated that BHLHE40 mRNA expression was significantly increased in MTB-1458-infected THP-1 macrophages in comparison to uninfected cells and was negatively correlated with the miR-378d expression trend ( Figure 6A). The interac-tion between miR-378d and BHLHE40 was identified by the dual-luciferase reporter using psiCHECK-BHLHE40-wild and psiCHECK-BHLHE40-mut plasmids ( Figure 6B). Results showed the Renilla/firefly ratio was significantly lower in miR-378d mimic than in negative control (NC) mimic in the psiCHECK2-BHLHE40-wild group, and the inhibition effect was dose-dependent ( Figure 6C), indicating that BHLHE40 is the downstream target of miR-378d. However, there was no interaction between miR-107 and BHLHE40 by the dual-luciferase reporter assay. The expression level of BHLHE40 and miR-378d was negatively correlated during MTB infection and BHLHE40 was a target of miR-378d. However, in-depth investigations are needed to elucidate the detailed regulation mechanisms.
107/miR-378d-BHLHE40 were presented as a loop in the networks, the function of these three miRNAs and their targets was further analyzed. A previous study reported that BHLHE40 was an essential repressor of interleukin (IL)-10 during MTB infection and may have a potential function of the anti-MTB mice model [13]. Therefore, this study chose the BHLHE40-miR-107/miR-378d-BHLHE40 regulation axis for further investigation.
qRT-PCR results indicated that BHLHE40 mRNA expression was significantly increased in MTB-1458-infected THP-1 macrophages in comparison to uninfected cells and was negatively correlated with the miR-378d expression trend ( Figure 6A). The interaction between miR-378d and BHLHE40 was identified by the dual-luciferase reporter using psiCHECK-BHLHE40-wild and psiCHECK-BHLHE40-mut plasmids ( Figure 6B). Results showed the Renilla/firefly ratio was significantly lower in miR-378d mimic than in negative control (NC) mimic in the psiCHECK2-BHLHE40-wild group, and the inhibition effect was dose-dependent ( Figure 6C), indicating that BHLHE40 is the downstream target of miR-378d. However, there was no interaction between miR-107 and BHLHE40 by the dual-luciferase reporter assay. The expression level of BHLHE40 and miR-378d was negatively correlated during MTB infection and BHLHE40 was a target of miR-378d. However, in-depth investigations are needed to elucidate the detailed regulation mechanisms. The psiCHECK2, psiCHECK2-BHLHE40-wild, and psiCHECK2-BHLHE40-mut plasmids were cotransfected with NC mimic or miR-378d mimic. The relative luciferase activity means Renilla to Figure 6. BHLHE40-miR-378d-BHLHE40 regulation axis verification. (A) THP-1 cells were infected with MTB-1458 at a MOI of 10 for 12 h. Cells were washed twice and a medium with gentamycin (100 µg/mL) was added. Thereafter, total RNA was collected at 6 and 24 hpi. BHLHE40 expression level was evaluated by qRT-PCR. (B) The binding site of miR-378d with BHLHE40 3 -UTR mRNA or mutant of BHLHE40 3 -UTR mRNA. (C) Dual-luciferase reporter assay was executed in 293T cells. The psiCHECK2, psiCHECK2-BHLHE40-wild, and psiCHECK2-BHLHE40-mut plasmids were cotransfected with NC mimic or miR-378d mimic. The relative luciferase activity means Renilla to firefly ratio (Rluc/Fluc) in miR-378d mimic or NC mimic. KD represents knockdown. "***" represents p < 0.001, "+" represents presence, "−" represents absence.

Discussion
The difference in phenotype induced by virulent and avirulent strains is an effective approach to explore TB pathogenesis. Previous studies showed that, in the virulent MTB strain H37Rv-infected macrophages, 36 genes related to immune response regulation, chemokine secretion, and leucocyte chemotaxis were upregulated, and 30 genes associated with amino acid biosynthetic and energy metabolism, connective tissue development, and extracellular matrix organization were downregulated, compared to avirulent straininfected macrophage [14]. Virulent M. bovis and BCG can also lead to different immune responses associated with TB pathogenesis [2]. This study compared the expression profile of virulent MTB-1458-infected macrophages or avirulent BCG-infected macrophages to uninfected cells. A unique transcriptional network, including TFs and mRNAs of MTBvirulent-strain-mediated DE-miRNAs, was identified to comprehend the pathogenesis of these two strains, which might be favorable to prevent the development of TB.

miRNA Expression Profile Comparison between MTB-1458-Infected Macrophages and BCG-Infected Macrophages
Many studies have highlighted the role of miRNA during MTB infection, including modulation of autophagy, apoptosis, or inflammation, to affect the susceptibility or survival of MTB [15][16][17][18]. In this study, miRNA expression profiles of MTB-1458-infected THP-1, BCG-infected THP-1, and uninfected THP-1 were analyzed by high-throughput sequencing. MTB-1458 induced a stronger miRNA expression profile variation in the entire level compared to BCG-infected macrophages. In all 18 unique DE-miRNAs, only miR-146a-5p (miR-146a), miR-378d, miR-874-3p (miR-874), and miR-140-3p have been reported to modu-late the essential biological processes during MTB infection [19][20][21][22][23]. Thus, it is worthwhile to thoroughly research the function of the remaining miRNAs, which would contribute to finding new mechanisms during the interaction between the host and bacteria.
The candidate targets for DE-miRNAs were analyzed, and they were mainly enriched in the Foxo signaling pathway (5 of 81), TNF signaling pathway (4 of 81), PI3K/Akt signaling pathway (5 of 81), Rap1 signaling pathway (4 of 81), and C-type lectin receptor signaling pathway (3 of 81). The first four pathways were mainly related to apoptosis, consistent with the result of GO BP. Previous studies have shown that, compared to the avirulent strain, inhibiting cell apoptosis seems to be an effective measure of virulent strain to escape the host's antibacterial action [30]. Furthermore, C-type lectin receptors were important receptors that could identify pathogen-associated molecular patterns (PAMPs) of MTB. For instance, Mincle recognizes trehalose 6,6 -dimycolate (TDM) to induce proinflammatory factor production [31,32], whereas mannose receptor (MR) is the receptor of ManLAM that facilitates survival of MTB [33,34]. Studying the specific functions of the C-type lectin receptor signaling pathway is conducive for controlling the development of TB by dominating these receptors.

Key DE-miRNAs
According to the TF-miRNA-mRNA Network miRNAs affect the cellular environment and signals by targeting to mRNAs. Meanwhile, their expressions are also regulated by TFs. Therefore, TF, miRNA, and mRNA can form multiple feed-forward and feed-backward loops, providing a novel idea for studying the key role of miRNAs in the disease pathogenesis. TF-miRNA-mRNA regulation has been applied to study cancer pathogenesis and livestock breeding [35][36][37]. This study sought hub miRNAs differentially expressed only in the MTB-1458-infected group and their TFs and mRNA by constructing the TF-miRNA-mRNA network.
The upregulated DE-miRNA network consisting of 5 miRNAs, 12 TFs, and 14 targets was generated ( Figure 4A). miR-146a has been extensively researched in TB and mediates the immune response by targeting IRAK-1, TRAF-6, and PTGS2 [20,38]. However, the mechanism of miR-146a upregulated expression is still unclear. These results uncovered that miR-146a may be activated by TF-CTNNB1. The function of the remaining four miRNAs in TB has not been studied yet. miR-335-3p had the highest node degree, and JUND-miR-335-3p-MYB/SH3KBP1 was regarded as a key axis in all interactions. JUND is a member of the AP-1 family of TFs and can regulate transcriptional programs of cellular differentiation, proliferation, and apoptosis [39]. SH3KBP1 is implicated in numerous cellular processes, including apoptosis, cell adhesion, and regulation of clathrin-dependent endocytosis. MYB can influence cell survival through the PI3K/Akt signaling pathway. Since apoptosis is one of the differential pathogenic mechanisms during BCG and MTB infection, we suppose miR-335-3p may participate in the cell apoptosis process.
The downregulated DE-miRNA network comprised of 13 miRNAs, 11 TFs, and 67 targets ( Figure 4B). Three miRNAs (miR-378d, miR-874, and miR-140-3p) have been verified that could play important roles during MTB infection. miR-3184-5p had the highest node degree, and MYB-miR-3184-5p-BMF/IRF1/PML was regarded as a key axis composed of the higher degree node. Studies demonstrated that miR-3184-5p is a tumor suppressor gene that can inhibit cancer cell proliferation [40,41]. In addition, BMF is a member of the Bcl-2 family that is related to apoptosis. IRF1 participates in the IFN-γ and IL-12 signaling pathways. PML is relevant to IFN-γ signaling pathways and SUMOylation. Thus, miR-3184-5p may modulate the pathways mentioned above during MTB infection. In addition, some genes are both target and TF; they, hereby, interacted with one miRNA, making it possible to form a feed-forward or feed-backward loop. For instance, BHLHE40, CEBPB, and IRF1 all act as vital mediators during pathogen-host interactions [13,42,43]. This network analysis provides directions for further experimental exploration.
The difference of virulence between BCG and MTB can lead to variant immune responses, including immunomodulatory genes, cytokines, and miRNAs. Recently, hostdirected therapy (HDT) is considered to be one of the effective ways of treating TB in the future, and the strategies of using miRNA, cytokines, and key immunomodulatory genes as the HDT are being broadly discussed [44][45][46][47]. This study had carried out in-depth data mining of the profiles of miRNAs in THP-1 cells after the virulent and avirulent strain infection and found the key TF/miRNA/mRNA regulation axes by systematic network analysis. The analysis results revealed that cytokines and key immunomodulatory genes related to TB possessed a regulatory relationship with the DE-miRNAs. Therefore, it is reasonable to infer that some of the DE-miRNA can be used for HDT against TB. However, their specific functions need to be validated in a series of subsequent experiments. Moreover, the displayed DE-miRNAs in this research are only based on the THP-1 cell; subsequent verification tests are suggested to be performed with primary macrophages or in vivo. This will help to provide more sound and reliable theoretical knowledge for the development of potential druggable targets for HDT against TB.
In conclusion, this was the first comprehensive TF-miRNA-mRNA network analysis study based on global transcriptome analysis of virulent MTB-infected THP-1 cells or avirulent BCG-infected THP-1 cells and uninfected THP-1 cells. Unique DE-miRNAs in MTB-1458-infected groups as well as their matched targets and TFs were selected to generate the TF-miRNA-mRNA network for excavating key TF, miRNA, mRNAs, and axes. The results provided an excellent suggestion for further understanding the specific immune response between MTB and the host.
The infection experiment was performed as described previously [50]. Briefly, by using phorbol myristate acetate (40 ng/mL, Sigma-Aldrich) to differentiate THP-1 cells for 12 h. The mid-log phase cultures (MTB-1458 and BCG) were pelleted and resuspended in Hanks' Balanced Salt Solution (HBSS, Gibco). Bacterial were dispersed into a single-cell suspension by using a 22-gage syringe. The bacterial concentrations (CFU/mL) were estimated by OD 600 values. The multiplicity of infection (MOI, BCG, or MTB-1458: THP-1 cells) was 10. After 12 h, cells were washed twice to remove extracellular bacteria, and then a medium with gentamycin (100 µg/mL) was added. Thereafter, total RNA was collected at 6 and 24 hpi. There were six groups, including MTB-1458-6 hpi, MTB-1458-24 hpi, BCG-6 hpi, BCG-24 hpi, uninfected-6 hpi, and uninfected-24 hpi. All experimental procedures using live MTB were performed at the Biosafety Level 3 Laboratory at Huazhong Agricultural University.

Total RNA Extraction and Integrity Analysis
Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with the instructions of the manufacturer. RNA purity and integrity were checked by agarose gel electrophoresis and quantified by determining the absorbance of A260 on SmartSpec (Bio-Rad, Hercules, CA, USA).

miRNA Profiling by High-Throughput Sequencing
A Balancer NGS Library Preparation Kit for small RNA/miRNA (GnomeGen, San Diego, CA, USA) was used for small RNA cDNA library preparation. In accordance with the instructions of the manufacturer, 3 µg of total RNA was ligated to 3 and 5 adaptors, sequentially, reverse transcribed to cDNA and amplified by PCR. Then, 10% natural polyacrylamide gel electrophoresis was used to separate the whole library, as well as cut and elute the bands corresponding to miRNA insertion. The purified small RNA libraries were quantified for cluster generation and 36 nt single-end sequencing analysis using Illumina GAIIx (Illumina, San Diego, CA, USA).

Bioinformatic Analysis
The potential targets of miRNAs were predicted by miRanda [51] and miRDB [52]. The potential TFs of miRNAs were predicted by the TransmiR version 2.0 database [53]. GO classification and KEGG pathway enrichment analysis of the final candidate targets of selected miRNAs were performed in DAVID 6.8 [54] and KONAS 3.0 [55] online database with default settings, respectively. Jvenn was adopted to create the Venn diagram [56]. The associations among miRNAs, mRNAs, and TFs were constructed and visualized through Cytoscape version 3.6.0.

miRNAs Screening Differentially Expressed Only in MTB-1458 Groups
As shown in the flowchart (Figure 2), DE-miRNAs from the BCG-infected groups vs. uninfected groups and DE-miRNAs from MTB-1458-infected groups vs. uninfected groups were selected first. Then, DE-miRNAs only in MTB-1458 groups were selected. Second, DE-miRNAs with a transcript per million (TPM) value of <10 were rejected. Finally, the union of DE-miRNAs upregulated at 6 or 24 hpi, and the intersection of DE-miRNAs downregulated at both 6 and 24 hpi, were further selected.

miRNA and mRNA Expression Validation by qRT-PCR
miRNA and mRNA expression levels were analyzed by qRT-PCR as described previously [50]. Briefly, by using the HiScript II Q Select RT SuperMix kit (Vazyme, Nanjing, China) to reverse transcribe RNA (1 µg) into cDNA. For miRNA, using a miRNA 1st-Strand cDNA Synthesis Kit (GeneCopoeia, Carlsbad, CA, USA) to synthesize cDNA. qRT-PCR were performed with AceQ qPCR SYBR Green Master Mix (Vazyme) in ABI ViiA 7 (Applied Biosystems, Carlsbad, CA, USA). Meanwhile, U6 and β-actin acted as the internal reference, respectively. The 2 −∆∆Ct method was adopted to analyze the genes relative expression. Table S2 displayed the related primer sequences.

Identification of the Relationship between BHLHE40 and miR-378d with Dual-Luciferase Reporter Assay
The dual-luciferase reporter assay was performed as described previously [50]. The pairing region with miR-378d in partial 3 -untranslated region (3 -UTR) of BHLHE40 was amplified using the primer pair F, 5 -cctcgagCCAAACCAAGGTCTGAGAAATG-3 and R, 5 -ttgcggccgcCTCTTCAAAAACAGGAATACATTCA-3 . Overlapping PCR was applied to generate the mutated 3 -UTR with the primer pair of F, 5 -TATTTTTGCCGGTCTGTACTTGTT-3 and R, 5 -AACAAGTACAGACCGGCAAAAATA-3 .

Statistical Analysis
Data were expressed as mean ± standard deviation. ANOVA was used to analyze the statistical significance of the data by GraphPad Prism (version 8.0). * p < 0.05, ** p < 0.01, and *** p < 0.001 were used to indicate statistically significant differences.